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ABSTRACT 

All axisymmetric self-similar equilibria of self-gravitating, rotating, isothermal systems are identified by 
solving the nonlinear Poisson equation analytically. There are two families of equilibria: (1) Cylindrically 
symmetric solutions in which the density varies with cylindrical radius as i?~", with < a < 2. (2) 
Axially symmetric solutions in which the density varies as f{0)/r'^, where r is the spherical radius and 
is the co-latitude. The singular isothermal sphere is a special case of the latter class with f{6) = 
constant. The axially symmetric equilibrium configurations form a two-parameter family of solutions 
and include equilibria which are surprisingly asymmetric with respect to the equatorial plane. The 
asymmetric equilibria are, however, not force- free at the singular points r = 0, oo, and their relevance to 
real systems is unclear. For each hydrodynamic equilibrium, we determine the phase-space distribution 
of the collisionless analog. 

Subject headings: galaxies: structure — galaxies: kinematics and dynamics — stars: kinematics — 
stars: formation — ISM: kinematics and dynamics — hydrodynamics 



1. INTRODUCTION 

Many astrophysical objects consist of equilibrium con- 
figurations in which self-gravity is resisted by pressure (ve- 
locity dispersion) and centrifugal forces. A particularly 
simple example, which is well-suited for practical analysis, 
is the case of an isothermal fluid with a linear pressure- 
density relation, p = pel, where p is the pressure, p is 
the density, and Cs is the sound speed which is taken to 
be a constant. Equilibrium configurations of iso thermal 
syste i TLS have been studied as ra odels of g alaxies ( Toomrt 



the other parameter, _B, controls the symmetry of the solu- 
tions with respect to the equatorial plane. For i? = 1, the 
solutions are up-down symmetric. These solutions have 
been previously obtained by Toomre (1982) and Hayashi 
et al. (1982) and include the singular isothermal sphere 
[A = 2) and the cold Mestel disk {A -^ oo) as limiting 
cases. 

The solutions with B ^ 1 are asymmetric with respect 
to the eq uatorial plane. This contradicts Lichtenstein's 



1982 ; Binney fc Tremaine 1987 ; see also Richstone 1980 



theorem (Lichtenstein 1933; Wavre 1932| ) which states 



Monet et al. 1981 for some spec i al cases) and newly f ormed 
stars (e.g., Hayashi et al. 1982; Kiguchi et al. 1987| ). 



Although some analytical solutions have been published 
previously, no systematic study of isothermal equilibria 
has been presented so far. In this paper we classify and 
derive analytically all possible self-similar axisymmetric 
equilibria of a self-gravitating isothermal system. We find 
that there are only two distinct classes of equilibria: 

1. Cylindrically symmetric equilibria, in which all 
quantities, such as density, potential, velocity, etc., 
depend on the cylindrical radius, R, only; 

2. Axially symmetric equilibria, in which the quantities 
are functions of both the spherical radius, r, and the 
co-latitude, 9. 

The cylindrically symmetric solutions are rather simple: 
the density varies as p oc i?~^^"+^\ and the azimuthal ve- 
locity as flR = v^ = v^qR~'^, where the allowed range of 
n is — 1 < n < 0. 

The axially symmetric solutions are more rich. First, the 
density always varies as p ex 1/r^ and the rotation curve is 
flat, v,p = Vipo = constant. Second, these equilibria form a 
two-parameter family of solutions. One of the parameters, 
A, determines the rotation velocity, A = 2 + i;,^„n/c?, and 
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that a barotropic, self-gravitating equilibrium must have 
a plane of symmetry perpendicular to the axis of rotation. 
The paradox is resolved by noting that the solutions with 
B y^ 1 are not force- free at two singular points, r = and 
r = oo, where Poisson's equation is ill-defined. External 
forces have to be applied to the matter at the singularities 
to hold the system in equilibrium. This is proved rigor- 
ously for the case without rotation, A = 2, and is likely to 
be correct also for rotating solutions. 

The paper is organized as follows. In §0, we present a 
general analysis of the problem and classify all possible 
self-similar equilibrium configurations of self-gravitating 
isothermal systems. We derive cylindrically symmetric so- 
lutions of the equation of hydrostatic equilibrium in §0. 
We identify a two-parameter family of axially symmetric 
solutions in ^ and describe the properties of these solu- 
tions in §@. In §0 we discuss the non-rotating case and in 
§0 the thin disk limit. In ^ we determine the steady-state 
distribution function of collisionless stellar systems, and 
we summarize the conclusions in 99. 

2. GENERAL CONSIDERATIONS AND SELF-SIMILAR 
SOLUTIONS 

We perform all calculations in spherical coordinates 
(r, 0, (p), but we also occasionally consider the cylindrical 
radius, R = rsinO. 
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A hydrostatic equilibrium configuration satisfies the mo- 
mentum equation with vanishing time derivative, and 
Poisson's equation. By definition, the f- and 0-components 
( "hat" denotes unit vectors) of the velocity vanish in equi- 
librium. Also, by the condition of axial symmetry, all 
derivatives with respect to the toroidal angle, y?, vanish. 
Thus, we need to consider only the f- and ^-components 
of the momentum equation and Poisson's equation for the 
potential: 
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is the axisymmetric Laplacian operator, <f> is the gravita- 
tional potential, v^ is the toroidal component of the ve- 
locity due to rotation which is, in general, a function of r 
and 6, and G is Newton's constant. The radial coordinate 
and the density are taken to be dimensionless throughout 
the paper. 

2.1. Equilibrium of an Isothermal Gas 

We begin with a general treatment of the problem, with- 
out any assumption of self-similarity. We take the gas to 
be isothermal, i.e., the pressure and the density are related 
as follows, 

P = <?sP. (2) 

where Cs = const, is the sound speed. Since an isothermal 
system is a special case of a barotropic system, p = p(/9), 
the following generic propert ies immediate ly follow from 
the Poincare-Wavre theorem ( Tassoul 1978| ) : (i) the angu- 
lar velocity is constant on cylinders centered on the axis of 
rotation, i.e., in cylindrical coordinates (i?, i^, z), the veloc- 



ity -v^ 



,{R) is independent of z; (ii) the effective grav- 



ity can be derived from an effective potential; and (iii) the 
effective gravity is normal to the isopycnic (i.e., constant 
density) surfaces. Tassoul (1978) has derived general rela- 
tions between p, (j), and v^p in a rotating barotropic system; 
we briefly re-derive some of these results for completeness. 
Eliminating v^ between Eqs. (tta,b), we arrive at the 
following differential equation: 
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(3) 



Any function of the argument ln(rsin0) is a solution of 
this differential equation. Absorbing the logarithm into 
the definition of the function, we write the solution as fol- 
lows 

4){r,e) = ~cl\\-ip{r,e) ^u{rs\x\9), (4) 

where u is an arbitrary function of the cylindrical radius 
R = rsaiO. The function u is related to the rotation ve- 
locity as follows. 



\{r,e)=vl{R) = R 



dujR) 
OR ■ 



(5) 



This relation is obtained upon substituting Eq. (0) into 
Eq. ^. 

Next, we use Poisson's equation to relate the density 
and gravitational potential. Substituting Eq. (0) into Eq. 
( p^ and using Eq. (||) and the identity 



V^g u(rsin6') 



^ R 
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we obtain the following equation for the density distribu 
tion, 
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Given the rotation curve, v^{R). the solution of Eq. (o), 
together with Eqs. (g),(|4|), and (^, completely determines 
the solution. To proceed further, we need to make some 
simplifying assumptions. 

2.2. Self- Similar Solutions 

We now look for self-similar solutions of Eq. (g) . Let us 
assume that v^ is described by a power-law in i? = rsmO. 
Then the density distribution is also a power-law in r, and 
we write 

f{0) _ v^ 

'''^ " i?" ' 



P = Po- 



(7) 



where f{9) is an unknown function to be calculated, po 
and Vipo are normalization constants, and a and n are pa- 
rameters. Substituting Eqs. (M) into Eq. (0), we obtain 



^h[0d9 ''""^de ^^^(^V"(rsin(?)2»+2 



= AnGpo 

(^) 
This equation can be satisfied only when the powers of r 

on the various terms match. There are two possible cases: 

1. Cylindrically symmetric solutions. In this case, the 
first term on the left-hand-side of Eq. (||) vanishes 
identically and 2n + 2 = a. 

2. Axially symmetric solutions. In this case, a = 2 and 
n = and the term proportional to v^o in Eq- (||) 
vanishes. These solutions have flat rotation curves. 

The above list exhausts all possible cases. Thus, there are 
only two families of self-similar, axisymmetric equilibrium 
configurations of self-gravitating isothermal systems with 
permanent rotation. 

We now describe the two families of solutions in detail. 
In §0 we discuss the cylindrical solutions, and in §§0-0 
discuss the axially symmetric solutions. 
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3. CYLINDRICAL SOLUTIONS 

In this case, the first term on the left-hand-side of Eq. 

) vanishes identically, so that 2n + 2 = a. This condition 
orces / to be f{6) = (sin6')~". The self-similar solution 
for p in terms of R = r sin 6 is easily obtained from Eq. 
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We notice that this solution is physical only for n < 0; oth- 
erwise, the density is negative. The gravitational potential 
is determined from Eq. (0). We have, 



Vo 



2ttG i?2"+2 



-u{R) + (f>Q, 



(9b) 



where 4>q is a constant that defines the zero-level of the 
potential and the function u{R) is obtained by integrating 
Eq. dq) for a given power-law rotation curve. 
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iin^O, 
if n = 0. 



The solution is well behaved everywhere except at the axial 
singularity,^ where a Dirac (5- function, i.e., a central mass 
"wire," may be located. This additional mass density pro- 
portional to S{R) may be calculated using Gauss' integral 
theorem for the flux of the gravitational field through a 
surface, 



Js 



nds = -AttGM. 



(10) 



where Ms is the mass enclosed within a volume bounded 
by a surface S and n is the unit normal outward from the 
volume. We choose the surface S to be an axial cylinder 
of radius AR and length L, centered at i? = 0. Since (p is 
a function of R, [cf., Eq. (9b)], the gradient \/(l> is radial 
and is equal to 



\/(f)^R 



(2n + 2)c^ 
Ai? 



'^vo 



Ai?2"+i 



(11) 



where R is the unit vector along the cylindrical radius. 
Calculating the enclosed mass via Gauss' theorem and tak- 
ing the limit Ai? — > (remember, n < 0), we obtain the 
following expression for the linear density of the "central 



wire. 
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if n < 0, 
if n = 0. 



(12) 



Clearly, the mass at the singularity becomes negative for 
|n| > 1. Thus, the self-similar, cylindrically symmetric 
solutions exist only for n in the range — 1 < n < 0. The 
magnitude of the rotation velocity v^o is arbitrary. It is 
worthwhile to note that for n = — 1, the axial 5- function 
singularity disappears {fig — 0) and the solution, Eqs. (^, 
is well behaved everywhere. This case corresponds to solid- 
body rotation, n{R) = v^p/R = constant. 

4. AXIALLY SYMMETRIC SOLUTIONS 
The equation for p [Eq. (o) with a ~ 2 and n = 0] reads 



1 d 



sin0^1n/(0) 



AttGpq 



I{0). (13) 



4.1. Regular Solution 
It is remarkable that the nonlinear, second-order dif- 



ferential equation (O) may be solved analytically and a 
general two-parameter family of self-similar solutions can 
be found in a closed and explicit form. All details of this 
computation are presented in Appendix H. The solution 
is 



p{T,e) = 



A^ 



Btan^(6'/2) 



27rG (rsin6')^ [l -t- S tan'4(6l/2)]^ ' 



(14a) 



where A> and B > are two free parameters. We have 
restricted A to be positive, because a change of sign of A 
is identical to the replacement: B -^ 1/B. The gravita- 
tional potential calculated from Eq. (0) may be written in 
the form 
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(14b) 



The denominator in the logarithm is the contribution due 



to rotation, u{rsm6) 



Vo 



ln(rsin6') + const., as follows 



from Eq. (M). The requirement of regularity of the solu- 
tion, i.e., the continuity of p and (p and their derivatives 
everywhere, constrains one of the free parameters, 



A = 2 + vlJc^s>2 



(15) 



(see §4.2 for more details). The other free parameter, B, 
remains unconstrained. Note that the density distribu- 
tion depends on the rotation velocity through Eq. dla), in 
agreement with the Poincae- Wavr e theorem. For the spe- 
cial case B = 1, the solution (14a) reduces to the solution 
obtained by Toomre (1982) and Hayashi et al. (1982). ^ 

4.2. Axial Singularity and a Generalized Solution 

The solution (jlj) is smooth and well-behaved in the do- 
main < 9 < TT. However, the differential operator in 
Eq. (O) is ill-defined at = 0. Therefore, a singular so- 
lution, proportional to the Dirac delta-function, S{9), is 
allowed on the axis. This additional mass density may be 
calculated using Gauss' theorem, Eq. ([lO|). As a surface of 
integration, we choose a segment of a cone of revolution 
centered on the axis consisting of the following pieces, 

Si^{r = ro: < 6* < A9} , 

S'2 = {ro < r < ro + Ar; 9 = A9} , 

S3^{r = ro + Ar; < 61 < A9}, 



As one can see, the velocity gradient term drops out. 
That is, the density distribution appears to be indepen- 
dent of the rate of rotation, in apparent contradiction to 
the Poincare-Wavre theorem. As we shall see from the 
careful analysis below, this is not the case. 

^The singularity of the solution at S = is due to the singularity of the differential operator of Eq. ( 
^The parameter A = 2n + 2 in Toomre (1982) and A = 27 in Hayashi et al. (1982). 



The unit normals to 
n2 ~ 9, and n^ — f, 



where rg, Ar, and A9 are constants. 

the three surface pieces are rii = — f 

respectively. The gradient of the potential for small angles 

reads 

Pel {(3-A)cl 



V0: 



r9 

on the axis. 







(16) 



where (3 = 2 + v'^q/c^. Performing the integration over 
the surface S = Si & S2 (S S3, taking the hmit A0 -+ 0, 
and noticing that the radial distance, Ar, at 6* = is the 
length along the axis, Arj^^g = L, we obtain the linear 
mass density of a singular density distribution on the axis. 



dMs cl I V, 



-A 



(17) 



The requirement of a non-singular density, /i^ = 0, gives 
back the relation (15). In the more general case, we see 
that A is related to the rotation velocity and the mass den- 
sity at the axis. The mass on the axis cannot be negative, 
hence we have the constraint 



A<2 



^.^0 



Icl 



(18) 



We note that the radial part of the spherical Laplacian, 
V^g, is ill-defined at r = 0, so that the solution may also 
contain a point mass, <5(r), at the origin. However, using 
Gauss' theorem for a sphere of radius Ar -^ enclosing 
the origin, we find no singular mass can be "hidden" at 
r = 0. 

In the rest of the paper we consider only the regular so- 



lutions with /Xg 
[equation ([l5|)] 



0. These solutions satisfy A = 2 + v ^jc 



5. PROPERTIES OF THE AXIALLY SYMMETRIC 
SOLUTIONS 

5.1. General Properties 

Figures^ — H show typical examples of the self-similar 
solution (M ) . Llontours of equal density for different val- 
ues of the parameters A and B are displayed. As is seen, 
the parameter A determines the overall shape of the mat- 
ter distribution. As A — > the configuration tends to the 
cylindrically symmetric limit, and an axial singularity with 
/is ^ is always present [see Eq. (p^)]. Figure |^ shows an 
example with A = 0.7. The density does not vanish at in- 
finity at 9 = 0, TT. Solutions without an axial singularity, 
fj,s = 0, exist for A > 2, cf., Eq. dl5). The non-rotating 
limit, A = 2, yields confocal ellipsoids/spheres, as shown 
in Fig. (see §0). For A > 2, toroidal configurations with 
rotation are obtained, as shown in Figs, raand 3. Interest- 
ingly, as A ^ 00, the profile flattens and tends to a thin 
disk configuration, as shown in Fig. for A = 20. 

The parameter B is responsible for up-down asymme- 
try. For B = 1, the solutions are symmetric with respect 
to the equatorial plane. The solutions with B > 1 are 
shifted (distorted) upwards, those with B < 1 are shifted 
downwards. Note that solutions with B and 1/B are iden- 
tical except that they arc turned upside-down with respect 
to each other. Note also that the solutions with B ^ 1 
viola te Lichtenstein's theorem (Lichtcnstein 1933; Wavre 
1932), which proves the existence of an equatorial plane 
of symmetry. The resolution of this paradox is discussed 
below (§|). 

5.2. The Finite System Limit 

We now determine some global properties, such as the 
total mass, M, gravitational energy, W, and angular mo- 
mentum, L, of the axially symmetric self-similar solutions 
as functions of A and B. These quantities diverge unless 



a cutoff is imposed at large radii. The most natural way 
to do this is to assume that the system is immersed in a 
medium with finite pressure, Pcxt, but vanishing density. 
As follows from Eq. (g), the cutoff surface in this case is 
the iso-density surface with pc = Poxt/Cg. The equation 
for the cutoff surface is 



rM - Vg^WFc, 



(19) 



where g = cl/{2nG) and 



0(9) = 



A^ Bta.ii^{9/2) 

iin^e [l + Btaii^{9/2)y 



is the angular part of the function, p{r,9). Notice that if 
Ps 7^ 0, M and W still diverge due to the infinite mass lo- 
cated on the axis. So, we set /i^ = 0, which means that we 
are restricted to the solutions with A> 2. The total mass, 
gravitational potential energy, and angular momentum are 
defined as follows. 



M 



pdV, W 



p (f> dV, L = / p V X r dV. 



(20) 

The binding energy is then E = —{W + K), where the 
kinetic energy \s K = Afw^g/2. The surfaces of constant 
density and constant potential do not coincide due to ro- 
tation, unless A = 2. For this reason, the potential at the 
cutoff surface is a function of position, (pc = 4>{rc{9),9). 
To get rid of insignificant constants in the expression for 
W , we redefine the gravitational potential such that the 
maximum potential at the surface is zero. This constrains 
the constant (pQ. The value of ^m, the extremum point 
of (f) along the rc{9), is found to satisfy the equation, 
tan'^(6l,„/2) = 1/B. Then, we have 

0o = c2lnp,-i;2„lny^-i;2„ln(A/2). (21) 

As shown in §pl the solutions with B ^ 1 are not force- free; 
they include an external gravitational potential gradient. 
We do not include the external potential in our definition 
of the gravitational energy, W. We now calculate all the 
quantities (note, L^ is the only nonvanishing component 
of the angular momentum) . The details of the calculation 
are given in Appendix M. The results are 
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iz = 



where g = cl/{2TrG), ip{x) = r'(a;)/r(a;) is the digamma 
function, and 7 = — V'(l) = 0.57721 ... is Euler's constant. 
The dependences of M, W and Lz as functions of A are 
plotted in Fig. 0. All these increase with increasing A, i.e., 
as the rotation velocity Vcpo increases. The dependence of 



these quantities as a function of B is trivial; all are pro- 
portional to (S^/"^ + B^^/"^) and have minima at i? = 1, 
the symmetric configuration. 

We emphasize that the truncation of a system at a finite 
p = Pc breaks the self-similarity of the solutions. Thus, the 
above equations are approximate and should be used with 
care. Nevertheless, Eqs. ( p2| ) suggest that both the specific 
angular momentum of a finite configuration, Lz/M, and 
the gravitational binding energy per unit mass, W/M, are 
functions only of A and independent of B. The kinetic en- 
ergy per unit mass, i^^o/^i is also a function only of A. Fi- 
nally, if we consider the specific entropy, s = S/M, where 
S oc — Jy plnpdV, we see that S differs from the grav- 
itational energy by additive and multiplicative constants 
only. Thus s is again independent of B. The question 
then is: if energy and entropy are indistinguishable and 
independent of B, what determines which value of B is 
selected by a finite isothermal system with a given value 
of L, /M? The answer, as we show in the next section, is 



the bc pundary condition, namely, the magnitude of extcr- 
nal forces acting on the mass. 

6. THE NONROTATING LIMIT, A = 2 

We first prove that nonrotating solutions, yl — 2, consist 
of nested confocal ellipsoids. The solution (14a) reads 
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Btan2(6'/2) 



ttG (r sin 61)' [l + Btan2(6'/2)]' 



(23) 



The equation of an iso-density contour line is obtained by 
setting p{r, 6) — pi — constant. Upon straightforward 
trigonometric simplifications, we obtain 



r{e) = 



2ciB 



7rGpi(l + B)2 
a(l-e2) 
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l-B 

l + B 
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1 + e cos 9 ' 



(24) 



This is the equation of an ellipse with eccentricity, e, and 
semi-major axis, a, with the origin located at one of the 
focii. The three-dimensional iso-density surfaces are the 
surfaces of revolution about the major axis. Note that the 
case B = 1 corresponds to the singular isothermal sphere. 
Next, we demonstrate that the solutions with B ^ 1 are 
not force- free. Let us consider an equidcnsity (which is 
also an equipotential) surface with density pi. The gravi- 
tational potential inside the surface produced by the "out- 
side" mass, p < pi, is linear in the vertical coordinate, z, 
as shown in Appendix 0: 



$ - -z ^2^^ V^ (i?i/2 - B-1/2) 
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(25) 



As one can see, $ ex ^/pi- On the other hand, the mass in- 
side the equipotential is M = vg^l'^ j Jpl B^/^ + B''^/'^) 



I.e., M oc 1/^/pi, as follows from Eq. (23). The total 
gravitational force exerted by the "outside''^ mass on the 
"inside" mass, F^ = — MV$, is, thus, independent of pi, 
i.e., independent of the the equidensity surface chosen: 
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(26) 



where e ~ ~{B— l)/(i?+ 1). Taking an equipotential sur- 
face infinitely close to the origin, one can see that there 
is a net finite gravitational pull, Fg, (upwards for B > 1 
and downwards for B < 1) acting on the matter at the 
singularity, r = 0. This force is produced by the gravity 
of all the mass of the system. For the system to remain 
in equilibrium, there must be an equal and opposite exter- 
nal force, — Fg, applied at r = 0. A similar consideration 
shows that there must also be an extra force acting at in- 
finity or at the last equidensity surface of a finite system. 
These external forces are uniquely related to the asymme- 
try parameter, B, and, hence, determine the structure of 
the equilibrium. When B — 1, the external forces vanish 
and such an equilibrium configuration is force-free. While 
this analysis is restricted to the non-rotating case, A — 2, 
similar results should hold for any A. We have verified 
this numerically. 

As we mentioned earlier, the asymmetric solutions ap- 
parently contradic t Lichtenstein's theorem ( Lichtenstein 



1933; Wavre 1932), according to which any rotating body 



for which the angular velocity does not depend on z (which 
is true for a self-gravitating isothermal system) always pos- 
sesses an equatorial plane of symmetry which is perpen- 
dicular to the axis of rotation. There are two assumptions 
(among others) used in the theorem which are violated by 
our B ^ 1 solutions: (i) no exernal potential is allowed and 
(ii) the system is assumed to be bounded and the density 
is nowhere infinite. As we have demonstrated, an asym- 
metric configuration with B ^ 1 exists only when external 
forces act on the system and pull it apart, in violation of 
(i). These forces are applied at the positions where con- 
dition (ii) is violated: either p or r is infinite, and the 
solution is ill-defined. 

7. THE THIN DISK LIMIT, A ^ oo 

The solutions with A > 2 have rotation. The systems 
become fiattened as the rotation increases and they tend 
to a thin disk as A ~ v'^o/c^ -^ oo. In these solutions 
the parameter B determines the "inclination angle." We 
define the inclination angle (or the opening angle of the 
"equatorial cone" ) , dm , as the angle at which the density 
is maximum for a fixed radial distance, r. In the laige-A 
limit, we can neglect the sin 6 in the denominator in Eq. 
(pla|). We then find 
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(27) 



Taking into account that ar c cot (x) = 7r/2 — arccot(l/x) 
and expanding the solution ( [l4a| ) about 6m, for small an- 
gles — 9m = '&^ arccot(i3^^), we obtain 
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(28) 

For large A, the sech function is very sharply peaked. 
Thus, the height-to-radius ratio of the disk may be esti- 
mated as 



f = A. 






B + l/B 



(29) 



In the symmetric case, B 

„2 



1, this reduces to the standard 
result, H/R ~ cJ/v'^q <C 1. As c^ — > the disk becomes 
the cold Mestel disk having a flat rotation profile. 

Toomre (1999, private communication) has studied the 
dynamics of a truncated cold thin "conical" disk with 
A ^ oo, B ^ 1. Figure ^ shows his results for the evolu- 
tion of a conical disk which is at i = has a self-similar 
form calculated in this paper. No external forces are ap- 
plied at either the inner edge or outer edge of the disk. 
Because of the unbalanced gravitational force acting on 
the innermost ring of the disk, the core region is acceler- 
ated upward. The system evolves as a function of time 
in a self-similar fashion, as shown by a sequence of curves 
in Figure o which corresponds to equally spaced times in 
logi. The initial conical disk is thus destroyed within a 
dynamical time, showing that a tree conical disk is highly 
unstable. On the other hand, if a downward external force 
of appropriate magnitude is applied on the innermost ring 
of the conical disk and an equal upward force is applied 
on the outermost ring, the system is fully stable. 

8. EQUILIBRIUM OF COLLISIONLESS STELLAR SYSTEMS 

There is a close analogy between fluid (gaseous) self- 
gravitatin g equilibrium objects and s tellar coUisionless sys- 
tems (cf., Binney & Tremaine 1987). The analogy arises 
because a fluid system is supported against gravity by 
pressure gradients, while a stellar system is supported 
by gradients in the stress tensor which in many respects 
are like the pressure except that they can be anisotropic. 
There is a unique connection between the density and ve- 
locity fields of fluid equilibrium configurations and the dis- 
tribution functions of stars in their coUisionless analogs. 
For any stellar distribution function, /, the profiles of den- 
sity, streaming velocity, v^, velocity dispersion, etc., are 



the density profile is insensitive to the direction of rotation, 
it yields, upon inversion, a distribution function which is 
even in C^, i-e., /+ = f{£,C^) + f{£,—C^). Thus, there 
are, in principle, infinitely many distributions which pro- 
duce identical density profiles and differ by an arbitrary, 
odd in £z, function, /_ = f{£,£^) — f{£,—£^), deter- 
mined by the sense of motion of individual stars. For real 
stellar systems, /_ can also be determined, given the av- 
erage rotation profile, v^. Henceforth, we focus on the 
symmetric part of the distribution function. We omit the 
subscript "-I-," since this should not cause any confusion. 
To determine /, we need to express the density profile 
as a function of the cylin drical radius and gravitational 
potential. From Eqs. ( |l4b| ) and (|l^), it follows that 



p(i?,, 



R^-^e^p[~i<f>~^o)/cl 



(31) 



We now observe that p{R, 0) is independent of B. Hence 
f{£,C^) derived from it is also independent of B.^ Since 
f{£,C^) is unique (up to /_), the distribution function 
proposed by Toomre (1982) as the starting point of his 
analysis for the specific case of i? = 1, is, in fact, valid for 
any value of B and reads. 



/(£,/:j = /o/:r^exp -£/c: 



(32) 



where /o is a normalization constant. We can now write 
the phase-space distribution function of a steady-state, 
self-gravitating system of stars, /(x, v), in an explicit form 
as follows, 



/(r,0,v) = /op(r,0)exp 



calculated uniquely by taking the moments of /. Simi- 
larly, given /5, v^p^ etc., an equilibrium /(x,v) may be de- 
termined. In general, there is no guarantee that the stellar 
distribution function so obtained will be physical, i.e., non- 
nega tive over the entire six-dim ensional phase space (see, 
e.g., Binney & Tremaine 1987 for more discussion). Re- 
markably, however, the distribution finction is guaranteed 
to be positive for isothermal systems: the structure of an 
equilibrium of a self-gravitating isothermal gas is identi- 
cal to the structure of a coUisionless system of stars. We 
now determine the distribution function of a stellar sys- 
tem which corresponds to the hydrostatic equilibria found 

in§| 

We consider here the simple case when the distribution 
function depends on two classical integrals of motion, the 
energy and axial component of the angular momentum.^ 
Lynden-Bell (1962), and subsequently Hunter (1975), de- 
veloped a method to derive distribution functions from the 
density and mean velocity profiles. We introduce the in- 
tegrals of energy and angular momentum of a particle as 
follows. 



2di 



„2 



In I 



_ _ (33) 

Here c^ = w^ — v^ is the velocity dispersion of stars and 
jfj^ol is their mean circular velocity. We again notice a 
remarkable fact: the parameter B defines the shape of the 
gravitational potential, but does not affect the velocity 
distribution of stars in this potential. 

9. CONCLUSION 

In this paper, we analytically calculated all possible 
self-similar, axisymmetric equilibrium states of a self- 
gravitating, isothermal gas with rotation. We showed 
that there are two distinct classes of hydrostatic equilib- 
ria; namely cylindrically symmetric equilibria and axially 
symmetric equilibria. The axially symmetric solutions are 
more physical, since the matter density vanishes at in- 
finity. Among the axially symmetric solutions, we found 
equilibrium states which are asymmetric with respect to 
the equatorial plane. Such states satisfy Poisson's equa- 
tion, but are not force-free at the singularities, r = 0, oo. It 
is the external forces that support the asymmetric shape. 
This example shows that self-similar solutions should be 
treated with caution and checked for possibly unphysical 
boundary conditions that may be "hidden" at their singu- 
lar points. 

^According to Jeans' theorem, the distribution function of a steady state, self-gravitating system may be presumed to be a function of at 
most three isolating integrals. Usually, tw o of thp integrals arp thp pnergy and axial angular momentum, while, in general, there is no simple 
analytical form for the third integral (see, Binney & Tremaine 1987 for discussion). 

^More precisely, it depends on B via 0(r, d) entering the definition of C. 



£ = {vl + vl + vl) /2 + (/)(r, 9) - 0o, £, = \v^\ r sin0, 

where Vi are the components of the particle velocity. Since 



Are the asymmetric configurations likely to be realized 
in nature? Since real systems are finite, we should trun- 
cate our solutions at some pc and throw away the "out- 
side" mass. For a truncated solution to be valid, the grav- 
itational potential of the discarded "outside" mass has to 
be replaced by a suitable external potential. As Eq. (25) 
shows, the potential must have a linear gradient in z (for 
a nonrotating solution) of a magnitude determined by pc 
and B. Such a potential may be produced, for instance, 
by an external object located on the z axis at a distance 
large compared to the size of the system. Second, to keep 
the system at rest, another force of equal magnitude and 
opposite sign must act on the central core. Since the force 
must act on the core alone, it must be of non-gravitational 
origin, which, clearly, is not easy to arrange. One could 
imagine, for example, that the core is ionized (by radia- 
tion of a central source, for instance), while the material 
outside the core is neutral. Then, if an external magnetic 
field threads the ionized core, one could imagine the field 
pinning the core, but having no influence on the neutral 
outer material. The distortion from equatorial symmetry 



would then be determined by the strength of the gravi- 
tational attraction of the external object and the coun- 
terbalancing force from the tension and curvature of the 
field lines. Although this construction is rather artificial, 
it demonstrates that asymmetric, non-force-free equilibria 
may, at least in principle, exist in nature. If they do, and 
if the systems are isothermal, they will resemble some of 
the -B 7^ 1 solutions derived in this paper. 
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We also thank George Rybicki for useful discussions and 
Kristen Menou for help with the translation from French of 
the proof of Lichtenstein's theorem given by Wavre (1932). 
This work was supported in part by NASA grant NAG 5- 
2837 and NSF grant AST 9820686. 



APPENDIX 



METHOD OF SOLUTION OF EQ. (|13|) 
Here we demonstrate how a general solution of Eq. ( |l3| ) can be found analytically. The equation we analyze reads 



^-i^G^^^^^^^(^))-^(^)' 



(Al) 



where k = AttGpq/c^. First of all, we get rid of the constant term on the left-hand-side. We observe that 



sm dO \ do 



-/?, 



(A2) 



where /3 is a constant. Thus, upon introducing a new function g{6) = f(0)/ sin 9, equation (Al) becomes 



suie^Une^\ng{e)]+ng{e)^Q. 



(A3) 



We now observe that sin 6* can be put into the denominator as follows, 80/ smO — 91n |tan(0/2)|. Let us redefine again 
the function g{0) and change the independent variable as follows. 



e = ln 



tan ■ 



w{0=\r,g{e)^\n[l^^^ 
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Then the differential equation is greatly simplified and reads 



— ^O + '^e-(«)=0. 



(A5) 



One can now reduce this equation to a first-order differential equation which can be integrated. Since Eq. ( A5) contains 
no explicit dependence on ^, we introduce a new function p which is a function of the old function w as follows. 



, , dwif) , d'^w(£) dpiw) dpiw) dwif) , .dp(w) 

Upon this substitution, the resulting differential equation, pp'^j = — kc™, is readily integrated for p{w) to yield 

dw{0 



d^ 



j{w) = ±VCi-2Ke«'(«). 



Note that the above substitutions are equivalent to Howard's transformation used by Toomre (1982). 
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Here Ci is a first constant of integration wliich we assume, to be specific at the moment, to be positive (see discussion 
below). Making use of another substitution. 



viO 



- p'^CO 



the resulting differential equation, v'^ = ±vy/Ci — 2kv, can be integrated again. We arrive at the expression: 



/cr 



In 



y/Ci ~ 2kv{o + vc; 



±i^~C2), 



(A8) 



(A9) 



where C2 is a second constant of integration. We now resolve this equation with respect to v{^). There are two solutions, 
namely 



where we used the short-hand notation e* — exp (±^/U[[£^ — C'2])- Comparing Eqs. (A4) and (A8), one can see that the 
density p{r, 9) ex J [9) ex u(^). Thus, for the den sity to be positive, the function u(^) must be real and positive^ also. This 
condition is satisfied for the upper sign in Eq. (AlO). Thus, we have 
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(All) 



Recalling all the definitions we made throughout the calculation, see Eqs. (A4) and (A8), we finally obtain the general 
solution of the differential equation (|l3| ) , 



f{e) = 



Ci/ {2Ksh?9) 



cosh^ (\/CT [In |tan(6l/2)| - C2] /2) ■ 



(A12) 



Clearly, since cosha; has no zeros for —00 < a; < 00, this solution is well behaved for all 6''s from the range: < < tt. 
Therefore, we investigate some properties of this solution in more details in the main text. Quite interestingly, solution 
(A12) can be simplified further and finally reads as follows. 



/W 



2^2 



B|tan(6l/2)|'' 



KSm 



1 + S|tan(6'/2)|' 



(A13) 



where A = \fC\ and B = exp (— VCaC2) > are new constants of integration. The "absolute value" signs may be 
omitted if < < tt, i.e., the tangent is non- negative. 

It can be rigorously shown that the particular choice of Ci > and C2 being a real number is unique, provided the 
solution to be physically meaningful. The derivation following Eq. (A7) holds for the complex- valued constants Ci and 
C2. Separating real and imaginary parts in Eq. (|A1C| ) and requiring u to be a real-valued function for all real ^, one can 
arrive, upon straightforward but cumbersome mathematical manipulations, at the conclusion that Ci must be positive 
and C2 may be complex, such that Im.y/CiC2 — Trfc, k — 0, 1,2, . . . and its real part is arbitrary. From this, it follows 
that e* > for even fc 's an d e* < for odd fc's. The additional constraint that v must be positive uniquely selects the 
solution given by Eq. ( All ). 

EXPRESSIONS FOR THE MASS, GRAVITATIONAL POTENTIAL ENERGY, AND ANGULAR MOMENTUM 
From Eqs. (14), using Eqs. (|l5| ) and (|l9|), we straightforwardly calculate 

27r(73/2 



M 



■I, 



W 



2 i"ff 



3/2 



In p, + 2 - ^ - (^ - 2) fin ^^, - 1 



I- --^]J 



T '^9 v 

Pc 



where g = c2/(27rG') and T, J', and /C are the integrals: 
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Changing the integration variable to z = Btan'^{9/2), the above integrals reduce to 
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The first integral, T, is symbolically identical to / = /„ z°'~^R{z)dz, where R{z) is a rational function. It can be 
evaluated in the complex plane. The branch cut is taken to be from to +oo along the real axis. We choose the contour 
of integration which consists of four pieces: (i) a circular part around the branching point z — of radius i?< — > 0, (ii) a 
linear part, i_|_, along the cut from above extending from _R< to i?>, (iii) a circular part of radius -R> -^ oo, and (iv) a 
linear part, L_, going along the cut from below and closing the contour. One can see that the integrals over the circular 
pieces vanish when appropriate limits are taken, the integral along L+ is the integral we arc looking for, and the integral 
along L_ differ from it by a phase shift. Then, it becomes 
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where res^j. denotes the residue of a function at a pole z^ and sum goes over all poles of function R{z). Applying this 
method, we obtain 



I = 



{A' - 4) 



16 cos(7r/A) 



(«' 



M+^-iM 



(B5) 



The second integral, J , consists of three terms, due to the logarithm of a product. The first of them is identical to /, 



I.e., Ji = Jo 



R{z) In A^ dz ^ I In A^ . The second one, J2 = /„ z" ^R{z) In z dz, may be evaluated either directly via 



residues, as above, because In z does not affect the convergence of the integral, or by noticing that J^ z" ^R{z) In z dz 
■^ /o°° z"~^R{z) dz and using the previous result. The evaluation of the last integral, J3 = /g°° z"~^(l + z)^^ ln(l + z) dz, 



is trickier because there is another branching point, z 
modify it as follows. 



-1, which coincide with the pole of the rational function. We 
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/3=3 



where B{m,n) = r(?Ti)r(n)/r(TO + n) is the bet a- function, ijj{x) = r'(a;)/r(x) is the digamma function, and r(a;) is the 
gamma- function. Using the identity 'i}){x + 1) = tpi^) + I/2; to express tp{S) and combining all terms together, we arrive 
at the following result, 
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where 7 = — V'(l) = 0.57721 ... is Eulcr's constant. This equation may be further simplified using the identities: ili{x+l) 
ip{x) + 1/x and ■(/; (^ + a;) — -f/) (i — a;) = 7rtan(7ra;). We arrive at the following result, 
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' ,0 4^2 TT /I 1 
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The last integral, /C, is similar to 2. It is evaluated via the residue theorem to yield 
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GRAVITATIONAL POTENTIAL OF A TRUNCATED, NONROTATING SYSTEM 

Here we analytically calculate the gravitational potential of the nonrotating, A — 2, system truncated along the cutoff 
equidensity (and, hence, equipotential) surface with density p = pc. Note, since we cut the system along the equipotential 
surface, the potential produced by the "inside" mass in the outer space and the potential inside the cavity due to all 
"outside" mass are of equal magnitude and add up to a constant. The direct evaluation of the gravitational potential via 
the volume integration involves elliptic integrals and turns out to be cumbersome. We use a different approach. 



Let us assume, without loss of generality, that B > 1. Then, the ellipsoids are shifted upwards, as follows from Eq. (23) 
and plotted in Fig. 0. We label each ellipsoid with the quantity 



A = ae = Ao/Vp, 



(CI) 



where a is the major axis, e is the eccentricity, and Aq — — a/cj/SttG (_B^/^ — i?^^/^). We now consider two such confocal 
ellipsoids; the quantities referred to a larger one are denoted by the "prime" , as shown in Fig. 0. Let us now fill the space 
between the two ellipsoids with matter of a homogeneous density p. The gravitational potential in the empty space inside 
the smaller ellipsoid is equal to the potential inside the large homogeneous ellipsoid, $|jjj(x), less the potential inside the 
small one, $int(x), having the sa me density. The gravitatio nal potential in the interior of a homogeneous prolate ellipsoid 
centered at the origin is known (Binney & Trcmaine 1987): 



$int(x) = -7rGp /62-^A,x? 
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where b = a\/l — e^ is the minor axis, x = (x, y, z), and 
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Taking into accomit that our coordinate system is shifted through A < 0, we write, ^ Aixf — AiR^ + ^3(2 + A)^. The 
potential inside the shell bounded by the surfaces A' and A is, thus, 
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Making A' infinitesimally close to A, we obtain the contribution to the potential due to the infinitely thin shell of 
"thickness" SA = A' - A and density p = (Ao/A)^: 



(^^insidc = -'2TTG 
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The gravitational potential inside the truncated system extending from the cutoff surface Ac — AqJJjTc to infinity is 
simply the integral over all thin shells with A > Ac, i.e., <I>cav = /<5*l^inside- The first term in Eq. ( |C5| ) evaluates to a 
constant, which may be absorbed into $0, the zero-level of the potential. The second term yields the coordinate dependent 
contribution we are looking for. 



<I'cav(x) = -z v/^G^ (i?i/2 - B-1/2) A3(e)V^, 



(C6) 

where e = (1 — i3)/(l + B) and ^3(6) = ^3(— e) is given by Eq. (|C3|). We emphasize that the potential is linear in the 
vertical coordinate, z, which corresponds to the homogeneous gravitational field. The strength of this field is determined 
by B and the cutoff density, pc, only. The gravitational acceleration, g = — V$cav, is upwards for B > \ and downwards 
for B <\. 
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Fig. 1. — Iso-density contours for A = 0.7, B = 1 and B = 3. The horizontal axis corresponds to R and the vertical axis to z. 



Fig. 2. — Same as in Fig. 
Fig. 3. — Same as in Fig. 
Fig. 4. — Same as in Fig. 



\ ior A = 2, B = 1 a.nd B = 3. 
^ior A = 3, B = 1 a.nd B = 3. 
]\ior A = 20, B = 1 and B = 10^ . 



Fig. 5. — The total mass M (solid curve), gravitational energy W (short-dashed curve), and angular momentum L^ (long-dashed curve), 
of a finite system vs. A. All quantities are in arbitrary units. 

Fig. 6. — The evolution of a free axially symmetric "conical disk" as calculated by Toomre (1999, private communication). Radial profiles 
of density in cylindrical coordinates are shown at times: t = (straight line) and (from botton to top) t = 0.25, 0.5, 1, 2, 4. The coordinates 
and time are in arbitrary units. 

Fig. 7. — Two confocal ellipsoids, used to calculate the potential of a nonrotating truncated system. 
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